clear
% x1=fft(1:5:320);
% x2=fft(2:5:320);
% x3=fft(3:5:320);
% x4=fft(4:5:320);
% x5=fft(5:5:320);
% y=fft(1:320);
% X=[x1;x2;x3;x4;x5];
% y_=complex(zeros(5,64));
% for i=1:64
% 
%     w=[];
%     for j=1:5
%         w=[w;1,cir((i-1+64*(j-1)),320),cir((i-1+64*(j-1))*2,320),cir((i-1+64*(j-1))*3,320),cir((i-1+64*(j-1))*4,320)];
%     end
%     y_(:,i)=w*X(:,i);
% end
% y_=reshape(y_.',1,[]);
% re=abs(y-y_)<0.0001;

% for i=1:15
%    a=(cir(i,64)*2^11); 
%    fprintf("%d, ",round(real(a)));
%     
% end














% w0=1;
% w1=cir(64,320);
% w2=cir(64*2,320);
% w3=cir(64*3,320);
% w4=cir(64*4,320);
% a1=x1(1)*w0+x2(1)*w1+x3(1)*w2+x4(1)*w3+x5(1)*w4
% a2=y(65)

% w0=cir(0,320);
% w1=cir(1,320);
% w2=cir(2,320);
% w3=cir(3,320);
% w4=cir(4,320);
% a1=x1(2)*w0+x2(2)*w1+x3(2)*w2+x4(2)*w3+x5(2)*w4
% a2=y(2)


% w0=cir(0,320);
% w1=cir(65,320);
% w2=cir(65*2,320);
% w3=cir(65*3,320);
% w4=cir(65*4,320);
% a1=x1(2)*w0+x2(2)*w1+x3(2)*w2+x4(2)*w3+x5(2)*w4
% a2=y(66)

% for i=1:80
%    a(i)=cir(i,320); 
%    fprintf("%.7e, ",imag(a(i)));
% end
% ar=real(a);
% ai=imag(a);

x0=1:320;
for i=1:5
x1((i-1)*64+1:i*64)=i:5:320;
end
x2=b2(x1,5);
x3=b2(x2,10);
x4=b2(x3,20);
x5=b2(x4,40);
x6=b2(x5,80);
x7=b2(x6,160);
% for i=1:320
%     fprintf("%d, ",x6(i)-1);
% end


S=fft(randn(1,320));
Srank=S(x7);
S1=complex(zeros(1,320));
S2=complex(zeros(1,320));
S3=complex(zeros(1,320));
S4=complex(zeros(1,320));
S5=complex(zeros(1,320));
S6=complex(zeros(1,320));
S7=complex(zeros(1,320));
for i=1:2:320
    [S1(i),S1(i+1)]=btf(Srank(i),Srank(i+1),cir(-0,2));
end
for i=1:4:320
    [S2(i),S2(i+2)]=btf(S1(i),S1(i+2),cir(-0,4));
    [S2(i+1),S2(i+3)]=btf(S1(i+1),S1(i+3),cir(-1,4));
end

for i=1:8:320
    for j=0:3
        [S3(i+j),S3(i+j+4)]=btf(S2(i+j),S2(i+j+4),cir(-j,8));
    end
end
for i=1:16:320
    for j=0:7
        [S4(i+j),S4(i+j+8)]=btf(S3(i+j),S3(i+j+8),cir(-j,16));
    end
end
for i=1:32:320
    for j=0:15
        [S5(i+j),S5(i+j+16)]=btf(S4(i+j),S4(i+j+16),cir(-j,32));
    end
end
for i=1:64:320
    for j=0:31
        [S6(i+j),S6(i+j+32)]=btf(S5(i+j),S5(i+j+32),cir(-j,64));
    end
end

for i=1:64
    [S7(i),S7(i+64),S7(i+64*2),S7(i+64*3),S7(i+64*4)]=btf5(S6(i),S6(i+64),S6(i+64*2),S6(i+64*3),S6(i+64*4),i-1);
end
out=real(S7)/320;
SS=[Srank;S1;S2;S3;S4;S5;S6;S7;];
pool=tablegen(S1);



S=fft(1:320);
S1=complex(zeros(1,320));
S2=complex(zeros(1,320));
S3=complex(zeros(1,320));
S4=complex(zeros(1,320));
S5=complex(zeros(1,320));
S6=complex(zeros(1,320));
S7=complex(zeros(1,320));

for i=1:160
    [S1(i),S1(i+160)]=btf_r(S(i),S(i+160),cir(i-1,320)); 
 end
for i=0:1
    for j=1:80
        [S2(i*160+j),S2(i*160+j+80)]=btf_r(S1(i*160+j),S1(i*160+j+80),cir(j-1,160)); 
    end
end

for i=0:3
    for j=1:40
        [S3(i*80+j),S3(i*80+j+40)]=btf_r(S2(i*80+j),S2(i*80+j+40),cir(j-1,80)); 
    end
end

for i=0:7
    for j=1:20
        [S4(i*40+j),S4(i*40+j+20)]=btf_r(S3(i*40+j),S3(i*40+j+20),cir(j-1,40)); 
    end
end

for i=0:15
    for j=1:10
        [S5(i*20+j),S5(i*20+j+10)]=btf_r(S4(i*20+j),S4(i*20+j+10),cir(j-1,20)); 
    end
end

for i=0:31
    for j=1:5
        [S6(i*10+j),S6(i*10+j+5)]=btf_r(S5(i*10+j),S5(i*10+j+5),cir(j-1,10)); 
    end
end


% w=[1,1,1,1,1;1,cir(-1,5),cir(-2,5),cir(-3,5),cir(-4,5);
%     1,cir(-2,5),cir(-4,5),cir(-6,5),cir(-8,5);
%     1,cir(-3,5),cir(-6,5),cir(-9,5),cir(-12,5);
%     1,cir(-4,5),cir(-8,5),cir(-12,5),cir(-16,5)];
w=[1,1,1,1,1;1,cir(1,5),cir(2,5),cir(3,5),cir(4,5);
    1,cir(2,5),cir(4,5),cir(6,5),cir(8,5);
    1,cir(3,5),cir(6,5),cir(9,5),cir(12,5);
    1,cir(4,5),cir(8,5),cir(12,5),cir(16,5)];
for i=0:63
    S7(i*5+1:i*5+5)=S6(i*5+1:i*5+5)*w.';
   %[S7(i*5+1),S7(i*5+2),S7(i*5+3),S7(i*5+4),S7(i*5+5)]=btf5(S6(i*5+1),S6(i*5+2),S6(i*5+3),S6(i*5+4),S6(i*5+5),indexr) 
    
end
S7=S7/320;












% s=1:8;
% s1=complex(zeros(1,8));
% s2=complex(zeros(1,8));
% s3=complex(zeros(1,8));
% 
% for i=1:4
%    s1(i)=s(i)+s(i+4);
%    s1(i+4)=(s(i)-s(i+4))*cir(i-1,8); 
% end
% for i=0:1
%    for j=0:1
%        s2(i*4+j+1)=s1(i*4+j+1)+s1(i*4+j+3);
%        s2(i*4+j+3)=(s1(i*4+j+1)-s1(i*4+j+3))*cir(j,4);
%    end
% end
% for i=1:4
%    s3(i*2-1)=s2(i*2-1)+s2(i*2);
%    s3(i*2)=s2(i*2-1)-s2(i*2);
% end
% s3


% S=fft(randn(1,16));
% Srank=S([1,9,5,13,3,11,7,15,2,10,6,14,4,12,8,16]);
% S1=complex(zeros(1,16));
% S2=complex(zeros(1,16));
% S3=complex(zeros(1,16));
% S4=complex(zeros(1,16));
% S5=complex(zeros(1,16));
% S6=complex(zeros(1,16));
% S7=complex(zeros(1,16));
% for i=1:2:16
%     [S1(i),S1(i+1)]=btf(Srank(i),Srank(i+1),cir(-0,2));
% end
% for i=1:4:16
%     [S2(i),S2(i+2)]=btf(S1(i),S1(i+2),cir(-0,4));
%     [S2(i+1),S2(i+3)]=btf(S1(i+1),S1(i+3),cir(-1,4));
% end
% 
% for i=1:8:16
%     for j=0:3
%         [S3(i+j),S3(i+j+4)]=btf(S2(i+j),S2(i+j+4),cir(-j,8));
%     end
% end
% for i=1:16:16
%      for j=0:7
%          [S4(i+j),S4(i+j+8)]=btf(S3(i+j),S3(i+j+8),cir(-j,16));
%      end
%  end
% 
% out=real(S4)/16;
% SS=[Srank;S1;S2;S3;S4];

pool=tablegen(S1);






function pool=tablegen(y)

pool(1).num=0;
pool(1).index=[];
pool(1).ri=[];
pool(1).zf=[];

for i=1:length(y)
    flag=0;
    for j=1:length(pool)
        if abs(real(y(i))-pool(j).num)<0.0001
            pool(j).index=[pool(j).index,i];
            pool(j).ri=[pool(j).ri,1];
            pool(j).zf=[pool(j).zf,1];
            flag=1;
            break;
        end
        if abs(real(-y(i))-pool(j).num)<0.000001
            pool(j).index=[pool(j).index,i];
            pool(j).ri=[pool(j).ri,1];
            pool(j).zf=[pool(j).zf,-1];
            flag=1;
            break;
        end
    end
    if flag==0
       tmp=length(pool)+1;
       pool(tmp).num = real(y(i));
       pool(tmp).index=[i];
       pool(tmp).ri=[1];
       pool(tmp).zf=[1];
    end
    
    flag=0;
    for j=1:length(pool)
        if abs(imag(y(i))-pool(j).num)<0.000001
            pool(j).index=[pool(j).index,i];
            pool(j).ri=[pool(j).ri,2];
            pool(j).zf=[pool(j).zf,1];
            flag=1;
            break;
        end
        if abs(imag(-y(i))-pool(j).num)<0.000001
            pool(j).index=[pool(j).index,i];
            pool(j).ri=[pool(j).ri,2];
            pool(j).zf=[pool(j).zf,-1];
            flag=1;
            break;
        end
    end
    if flag==0
       tmp=length(pool)+1;
       pool(tmp).num = imag(y(i));
       pool(tmp).index=[i];
       pool(tmp).ri=[1];
       pool(tmp).zf=[1];
    end
end

end



function [c,d]=btf(a,b,w)
    c=a+b*w;
    d=a-b*w;
end

function [c,d]=btf_r(a,b,w)
    c=a+b;
    d=(a-b)*w;
end

function [a,b,c,d,e]=btf5(f,g,h,y,j,indexr)
w=[];
for i=0:4
    index=indexr+i*64;
    w=[w;1,cir(-index,320),cir(-index*2,320),cir(-index*3,320),cir(-index*4,320)];
end
o=w*[f;g;h;y;j];
a=o(1);
b=o(2);
c=o(3);
d=o(4);
e=o(5);
end

function [a,b,c,d,e]=btf5_i(f,g,h,y,j,indexr)
w=[];
for i=0:4
    index=indexr+i*64;
    w=[w;1,cir(index,5),cir(index*2,320),cir(index*3,320),cir(index*4,320)];
end
o=w*[f;g;h;y;j];
a=o(1);
b=o(2);
c=o(3);
d=o(4);
e=o(5);
end



function x=b2(y,blk)
    blknum=length(y)/blk;
    for i=1:blk
       ypiece= y((i-1)*blknum+1:blknum*i);
       yout(1:blknum/2)=ypiece(1:2:end);
       yout(blknum/2+1:blknum)=ypiece(2:2:end);
       x((i-1)*blknum+1:blknum*i)=yout;
    end
end






%X1=fft([1,4]);
%X2=fft([2,5]);
%X3=fft([3,6]);

%[a,b,c]=buttf3(X1(1),X2(1),X3(1),6,0)
%[d,e,f]=buttf3(X1(2),X2(2),X3(2),6,1)



function [p,q]=buttf2(x,y,w)
    p=x+w*y;
    q=x-w*y;
end
function [p,q,r]=buttf3(x,y,z,N,M)
    n=N/3;
    mat=[1,cir(M,N),cir(M*2,N);1,cir(M+n,N),cir((M+n)*2,N);1,cir(M+n*2,N),cir((M+n*2)*2,N)];
    %mat=[1,1,1;1,cir(1,N),cir(2,N);1,cir(2,N),cir(4,N)];
    in=[x;y;z];
    out=mat*in;
    p=out(1);
    q=out(2);
    r=out(3);
end

function w=cir(k,N)
    w=exp(-1j*2*pi*k/N);
end

